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Abstract: 

Generalized linear models (GLMs) have been used widely for modeling the 
mean response both for discrete and continuous random variables with an 
emphasis on categorical response. Recently Yang, Mandal and Majumdar 
(2012b) considered full factorial and fractional factorial locally D-optimal 
designs for binary response and two-level experimental factors. In this 
paper, we extend their results to a general setup with response belonging 
to a single-parameter exponential family and for multi-level predictors. 
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1 Introduction 



Binary responses and count data are usually modelled using generalized linear models 
(GLMs). GLMs have been widely used for modelling the mean response both for 
discrete and continuous random variables with an emphasis on categorical response. 
Although the methods of analyzing data using these models have been discussed in 
depth in the literature (McCullagh and Nelder (1989), Dobson and Barnett (2008)), 
only a limited number of results are available for the choice of optimal design of 
experiments under GLMs (Khuri et al. (2006), Woods et al. (2006), Stufken and Yang 
(2012)). Recently Yang, Mandal and Majumdar (2012a) obtained some theoretical 
results on locally D-optimal factorial designs for 2 2 experiments with binary response. 
In a subsequent paper, Yang, Mandal and Majumdar (2012b) considered D-optimal 
factorial designs with k two-level predictors and binary response. In this paper, we 
extend their results to a much more general setup with multi-level predictors and a 
response belonging to a single-parameter exponential family. 
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Following Yang, Mandal and Majumdar (2012a, 2012b), we consider locally optimal 
designs, as Fisher's information matrix contains the unknown parameters, and so 
does the D-optimality criterion. The optimality criterion can be written in terms of 
the variances, or measures of information, at each of the design points. Note that 
these variances depend on the parameters through the link function. The locally 
D-optimal designs can be quite different from the uniform design which is commonly 
used in practice, especially when some of the regression parameters stay significantly 
away from zero (see also Yang, Mandal and Majumdar (2012b)). For illustration 
suppose we consider a linear predictor involving two two- level factors, namely, family 
history (xi) and physical exercise (22), which leads to the linear predictor 77 = 
/3q + [3iXi + (32X2- If we have solid reasons to believe that family history is a 
much more dominant factor than physical exercise, and on the basis of prior 
knowledge, an initial choice of parameters (/3q, {3±, $2) = (2,2,0.05) is reasonable, 
then the optimal allocation of 100 experimental units will be (6, 28, 33, 33) on the 
four distinct combinations of factor levels, namely, the design points. On the other 
hand, if the practitioner decides to use the popular uniform design (25, 25, 25, 25) 
with equal number of replicates per design point, then the relative loss of efficiency 
is about 5%. 

This paper is organized as follows. In section 2, we provide the preliminary setup for 
the generalized linear models under consideration. We discuss the characterization of 
locally D-optimal designs in section 3. In section 4 we discuss our search algorithms, 
both theoretically and numerically, for obtaining optimal approximate and exact de- 
signs. In section 5, we illustrate our results with two real examples. We conclude 
with some remarks in section 6. 



2 Preliminary Setup: Generalized Linear Models 



In this paper, we consider experiments with a univariate response Y which follows an 
exponential family with a single parameter 9 in canonical form 

f(y;9) = ex V {yb(9) + c(9) + d(y)} , 

where b(9),c(9) and d(y) are known functions. For examples, 

binary Y ~ Bernoulli (9) 
count response Y ~ Poisson(0) 

Y ~ Gamma(fc, 9), k > fixed 
y ~ N(9,a 2 ),a 2 > fixed. 



' exp {y\og-^ + log(l - 9)} , 
exp {y log 9-9- logy!}, 

I-Hogfl + logf^} 



exp < y- 



exp{^-^-^-ilog(2 7 ra 2 )} 



Consider independent random variables Y±, . . . , Y n with corresponding covariates or 
predictors Xi,...,x„, where x, = (xa, . . . , Xid)' G Under a generalized linear 
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model, there exists a link function g and parameters of interest (3 = (/3i, . . . , /3d)', 
such that 

E(Yi) = fii and 77* = #(/^) = xf/3. 

It is known that, for example, see McCullagh and Nelder (1989) or Dobson and 
Barnett (2008), 



var(Fj) 



provided b(8) and c{9) are twice differentiable. Since Y^s are independent, the d x d 
Fisher information matrix I = (Ijk) can be obtained by 

/ <9l dl_ \ _ X^fc / fytj 
Jfc ~ 5/3 J ~ ^ var(F i ) ^ % 

where Z = log f(Y; 9) could be written as a function of (3 too. 

Suppose there are only m distinct predictor combinations x 1; . . . , x m with numbers 
of replicates ni, . . . , n m , respectively Then 

ife ~ ^ var(^) ^77 J - n ^ Y ^{Y i )\dr ]i ) %lk 
where Pi = n^/n, % — 1, . . . , m. That is, the information matrix I can written as 

I = nX T WX 

where X = (xi, . . . , x m ) T is an m x d matrix, and W = diag(piu>i, . . . ,p m w m ) with 

1 rdfnY 

var(Fj) \drji J 
For typical applications, I or X T WX is nonsingular. 

Suppose the link function g is one-to-one and differentiable. Further assume that ^ 
itself determines var(Yj), that is, there exists a function h such that var(Yj) = h(rji). 
Let v = ((g" 1 )') /h. Then Wi = u(rji) = v (x//3) for each z. We illustrate it using 
the four cases below. 

Example 2.1: Binary response Suppose Yi ~Bernoulli( / u.j). Then h = g~ 1 (l—g~ 1 ). 
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For commonly used link function g, 
1 



1/(77) 



*(»7)[l-*(i7)] 



(expje"} - 1) [log (1 - exp{-e"})] 5 



exp{2rj— e v } 
k 1— cxp{— e 7 !} ' 



for logit link; 
for probit link; 

for complementary log-log link; 
for log-log link. 



Example 2.2: Poisson count Suppose Y> 
function g = log, h = v = exp and Wi = v{j]i) 



j Poisson(/ij). For the canonical link 
exp{7/j} for each i. 



Example 2.3: Gamma response Suppose Y{ ~ Gamma(fc, Hi/k) with known k > 
0. For the canonical link function g(p) = 1/fi, h(rj) = l/(krj 2 ), i>(rj) = k/rj 2 and 
Wi = v(r)i) = k/rj 2 . A special case is k — 1 which corresponds to the exponential 
distribution. 



Example 2.4: Normal response Suppose Yi 
the canonical link function g(fi) = fi, h{rf) = a 2 



~ N(fii, a 2 ) with known a 2 > 0. For 
1/(77) = 1/cr 2 and Wi = 1 /a 2 . 



3 Locally D-optimal Designs 



In this paper, we consider experiments with a design matrix X = (xi, . . . , x m ) r which 
consists of m distinct experimental settings, that is, m distinct combinations of values 
of d covariates. The response Yi follows a single-parameter exponential family under 
the generalized linear model 77^ = g(E(Yi)) = x.J/3. For instance, in a 2 3 experiment 
with binary response, 77 = (3 + (3ix± + /3 2 x 2 + /?3X 3 + (323X2X3 represents a model that 
includes all the main effects and the two-factor interaction of factors 2 and 3. The 
aim of the experiment is to obtain inferences about the parameter vector of factor 
effects (3; in the preceding example, (3 = (/3o, (3i, 02, 03, ^23) • 

In the framework of locally optimal designs, we assume that Wi = z/(x^/3), i = 
l,...,m, determined by regression coefficients (3 and link function g, are known. 
For typical applications, Wi > 0, i — 1, . . . , m. Then the design problem we consider 



here is to find the "optimal" allocation of n experiment units into (n±, 
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such 



that n± + • • ■ + n m = n, known as exact design, or alternatively, find the "optimal" 
proportion pi = rii/n, i = l,...,m, known as approximate design. In this paper, 
we mainly focus on the approximate designs, that is, p = (pi, ■ ■ ■ ,p m )' such that 
Px + • • • + Pm = 1- Since the maximum likelihood estimator of (3 has an asymptotic 
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covariance matrix that is the inverse of nX'WX, a (locally) D-optimal design is a p 
which maximizes 

/(p) = \X'WX\. 

For example, for the 2 2 experiment with main-effects model and binary response, 
/(p) = \%(p 1 p 2 pzwiw 2 w^ + p 1 p 2 piWiw 2 w A + piPzPiWxWzWi + p 2 p3P4W 2 w 3 w 4 ). See 
Yang, Mandal and Majumdar (2012a) for a detailed analysis of this problem. For 
general cases, the following lemma expresses the objective function /(p) = |XWX| 
is an order-<i homogeneous polynomial of pi, . . . ,p m . This is useful in determining 
the optimal p^s. See Gonzalez-Davila, Dorta-Guerra and Ginebra (2007) and Yang, 
Mandal and Majumdar (2012b). 

Lemma 3.1 Let X[ii,i 2 , . . . ,id] be the d x d sub-matrix consisting of the i\th, . . ., 
ijth rows of the design matrix X . Then 

f{v) = \X'WX\= \X[h,...,i d ]\ 2 ■Pi 1 w il ---p id w id . 

l<il<-<i d <m 

Based on Lemma 13. 1\ Yang, Mandal and Majumdar (2012b) developed characteri- 
zation theorems for D-optimal designs and saturated D-optimal designs for two-level 
factors with Xy G {—1, 1} (Yang, Mandal and Majumdar (2012b), section 3.1). Note 
that the entries of X in Lemma 13.11 can be any real number, which allows us to ex- 
tend their results to multiple-level factors. Following Yang, Mandal and Majumdar 
(2012b) we define for each i — 1, . . . , m, 

f I — z 1 — z 1 — z 1 — z \ 

/*(*)=/( i Pi, Pi-l,Z,~ Pi+i,...,- p m ), 0<Z<1. (1) 

V 1 - Pi 1 - Pi 1 - Pi 1 - Pi / 

Note that fi(z) is well defined for all p satisfying /(p) > 0. Lemma 7.1 and Lemma 7.2 
in Yang, Mandal and Majumdar (2012b) could be applied for our case too. Thus we 
obtain theorems below to characterize locally D-optimal allocations and saturated 
D-optimal allocations for experiments with multiple-level factors and fairly general 
response under GLMs. Note that the design matrix X here consists of m x d real- 
number entries, which generalizes the results in Yang, Mandal and Majumdar (2012b). 

Theorem 3.1 Suppose f (p) > 0. Then p is D-optimal if and only if for each 
% — 1, . . . , m, one of the two conditions below is satisfied: 

(i) p i = 0andf i (l) < ^/(p); 
(n) 0<p l <±andf l (0) = 1 ^f(p). 

Theorem 3.2 Let I = {z'x, . . . , id} C {1, . . . , m} be an index set satisfying \X[ii, . . . , 
id}\ 7^ 0. Then the saturated design satisfying p ix = p i2 = ■ ■ ■ = p id = h is D-optimal 
if and only if for each i 1, 



E 



W{Qui\Q}]| 2 < \x[i u i 2 ,...,, d }\- 



j el 3 1 
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Theorem 13.11 and Theorem 13.21 are for X with real-number entries and thus can be 
applied with multiple-level factors. We illustrate the results using the two examples 
below. 



Example 3.1: One two-level covariate and one three-level covariate: Suppose 
the design matrix X is given by 



X = 



( 1 


1 


1 \ 


1 


1 





1 


1 


-1 


1 


-1 


1 


1 


-1 





^ 1 


-1 


- 1 ) 



Then, for example, p\ = p 2 = Pa = 1/3 is D-optimal if and only if v\ + Av 2 < v 3 , 
V\ + v 2 + v± < v 5 and 4vi + 4t> 2 + v 4 < v 6 , where Vi = l/wi, i = 1, . . . , 6. On the 
other hand, p 3 = p 5 = p 6 = 1/3 is D-optimal if and only if v 3 + 4v 5 + 4t> 6 < V\, 
v 3 + v 5 + v e < v 2 and Av 5 + v e < V4. 



Example 3.2: 2x3 factorial design: Suppose the design matrix 



X 



( 1 
1 
1 
1 
1 

\ 1 



1 

-1 
1 


-1 



1 \ 

-2 
1 
1 

-2 
1 / 



(2) 



where the four columns correspond to effects J, A (two- level), Bi (three-level, linear 
component), and B q (three-level, quadratic component). Then, for example, p\ = 
P2 = P3 = P4 = 1/4 is D-optimal if and only if V\ + v 2 + V4 < v 5 and V\ + V3 + v 4 < 
vq] P2 = P3 = P4 = P5 = 1/4 is D-optimal if and only if v 2 + V4 + v$ < V\ and 
v 2 + ^3 + v 5 < Vq; ps = Pi = p$ = pe = 1/4 is D-optimal if and only if v 3 + f 4 + v 6 < v 1 
and v 3 + v 5 + v % < v 2 . 



4 Searching for D-optimal Designs 



Yang, Mandal and Majumdar (2012b) developed very efficient algorithms searching 
locally D-optimal approximate designs or exact designs with binary response and two- 
level factors. Essentially the same algorithms could be used for maximizing /(p) = 
|XTVX| for more general setup, as long as X consists of m distinct <i-dimensional real- 
number vectors and the response belongs to a single-parameter exponential family. 
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4.1 Lift-one algorithm for D-optimal approximate design 

For finding p maximizing /(p) = |XWX|, the lift-one algorithm, as well as Lemma 7.1 
and Lemma 7.2 in Yang, Mandal and Majumdar (2012b), could be applied for our 
case, with some modifications. Recall that fi(z) is defined in equation ([1]). 

Lemma 4.1 (Yang, Mandal and Majumdar, 2012b, Lemma 7.1) Suppose f (p) > 0. 
Then for i — 1, . . . , m, 

f i (z) = az(l-z) d ~ 1 + b(l-z) d , (3) 

for some constants a and b. If pi > 0, b = fi(0), a = ~q~^t / otherwise, 
b = f ( p ), a = fi (±) • 2 d - b. Note that a > 0, b > 0, and a + b > 0. 

Lemma 4.2 (Yang, Mandal and Majumdar, 2012b, Lemma 7.2) Let l(z) = az(l — 
zf~ x + 6(1 - z) d with < z < 1 and a > 0,b > 0,a + b > 0. If a > bd, then 
max z l(z) = (j^Y' 1 at z = < 1- Otherwise, maxj(z) = b at z = 0. 

Lift-one algorithm (for multiple- level factors and single-parameter exponential fam- 
ily response) 

1° Start with arbitrary p = (pi, ■ ■ ■ ,p m )' satisfying < Pi < 1, i = l,...,m. 
Compute / (po)- 

2° Set up a random order of i going through {1,2,..., m}. 

3° For each i, determine fi{z) as in expression ([3]). In this step, either fi(0) or 
fi (|) needs to be calculated according to Lemma [4. 1[ 

4° Define pl° = (jE^Pi, • • • , izJrPi-i, z*, jE^Pi+i, ■ ■ ■ , ir^Pm) , wnere z * maxi- 
mizes fi(z) with < z < 1 (see Lemma H~2j) . Note that /(p* ) = fi(z*). 

5° Replace p with pi'\ / (p ) with /(p* ). 

6° Repeat 2° ~ 5° until convergence, that is, /(po) = /(p* ) for each i — 1, . . . , m. 

Convergence and performance of lift-one algorithm: To guarantee the con- 
vergence, we may modify the lift-one algorithm as in Yang, Mandal and Majumdar 
(2012b, Section 3.3.1). A similar proof could be applied here to show that (1) if the 
lift-one algorithm or its modified version converges at p*, then p* is D-optimal; (2) the 
modified lift-one algorithm is guaranteed to converge (Yang, Mandal and Majumdar, 
2012b, Theorem 3.3). 
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Yang, Mandal and Majumdar (2012b, Section 3.3.1) also compared the time cost of 
lift-one algorithm with commonly used nonlinear optimization algorithms including 
Nelder-Mead, quasi-Newton, conjugate gradient, and simulated annealing. Overall, 
Lift-one algorithm could be 100 times faster than those algorithms. In this paper, we 
run more simulations to check how the computational time and number of non-zero 
Pi's vary across the number of factors and range of parameters. 

Figure 1: Performance of lift-one algorithm for m = 2 k and d — k + 1 




Figure 1 shows the time cost and number of nonzero p^s on average based on 1000 
simulated (3 for 2 k main-effects model with logit link. The time cost is based on a 
Windows Vista PC with Intel Core2 Duo CPU at 2.27GHz and 2GB memory. The 
relationship between the time cost and the number of support points m = 2 k is close 
to linear for moderate k. The time cost on average is only 1.25 sees, 1.74 sees, and 
1.84 sees on average for m = 2 7 = 128 with /3j's follow iid U{— 3,3), U(— 1,1), and 
U(— 0.5, 0.5), respectively. In the meantime, the corresponding number of nonzero 
Pi's on average are 28, 48, and 67 respectively. 

From a practitioners point of view, it is often desirable to keep the number of support 
points of a design small. Due to Lemma 14. 2\ the lift-one algorithm may force a Pi 
to be exactly zero. It is an advantage over commonly used nonlinear optimization 
algorithms which may keep an unimportant pi a tiny value but never let it to be zero. 
In many cases, it is hard to distinguish "negligible" from "tiny". As the range of /3j's 
increases from U(— 0.5, 0.5) to U(—3, 3), the mean number of nonzero pj's reduces by 
one half due to more and more small lOj's are generated. 
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4.2 Exchange algorithm for D-optimal exact designs 

Given the total number of experimental units n, to find D-optimal integer allocation 
n = (rii, . . . , n m )', Yang, Mandal and Majumdar (2012b) proposed another algorithm 
for two-level factors and binary response, called exchange algorithm, which adjusts 
rii and nj simultaneously for randomly chosen index pair The essentially same 

algorithm after some modifications could be used for a general design matrix X 
consisting of real numbers. The goal is to find the optimal n which maximizes /(n) = 
|X'W n X|, where W n = diag{niu>i, . . . , n m w m }. To do this, we need a modified version 
of Lemma 7.5 and Lemma 7.4 in Yang, Mandal and Majumdar (2012b). 

Lemma 4.3 (Yang, Mandal and Majumdar, 2012b, Lemma 7.5) Fixing 1 < i < j < 
m, let 

fij( z ) = f(n 1 ,...,n i - 1 ,z,n i+ i,...,nj-.i,s-z,n j+1 ,...,n m ) 

= Az(s-z)+Bz + C(s-z) + D, (4) 

where s = n { + rij. Then (i) D > =>• B > and C > 0; (ii) B > or C > 
=>• A > 0; (Hi) f(n) > A > 0; (iv) D = f(m, . . . , n^i, 0, n i+1 , . . . , 
n j - 1 ,0,n j+1 , n m ). (v) Suppose s > 0, then A = ^ (2 fa (f) - fa(0) - fa(s)), 
B = \{fa{s)-D),C= 1(^(0) -D). 

Lemma 4.4 (Yang, Mandal and Majumdar, 2012b, Lemma 7.4) Let g(z) = Az(s — 
z) + Bz + C(s — z) + D for real numbers A > 0, B > 0, C > 0, D > 0, and integers 
s > 0, < z < s. Let A be the integer closest to sA+ 2 ^~ C ■ 

(%) J/0 < A < s, thenm&x < z < s g(z) = sC + D + (sA + B - C)A - AA 2 at z = A. 
(ii) If A < 0, then maxo< z < s g(z) = sC + D at z = 0. 
(Hi) If A > s, then maxo< 2 < s g(z) = sB + D at z = s. 

Exchange algorithm for D-optimal integer-valued allocations: 

1° Start with initial design n = (rii, . . . ,n m )' such that /(n) = > 0. 

2° Set up a random order of going through all pairs 

{(l,2),(l,3),...,(l,m),(2,3),...,(m-l,m)} 

3° For each let s = rij +rij. If s = 0, let n*- = n. Otherwise, calculate fa(z) 

as given in expression (j4j). Then let 

n ij = • • • ? z *i n i+l, • • • j n j-li S — Z*, . . . , n m ) 

where the integer z* maximizes fa(z) with < z < s according to Lemma [4.41 
Note that /(n£) = /^(^) > /(n) > 0. 

4° Repeat 2° ~ 3° until convergence (no more increase in terms of /(n) by any 
pairwise adjustment). 
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5 Real Examples 



We use two real examples to show how our results work. 

Example 5.1: Paint gun: The following example is taken from Wu and Hamada 
(2009) and has been modified to suit our purpose. Hale-Bennett and Lin (1997) 
reported an experiment that was performed to improve a painting process of charcoal 
grill parts, which are painted with paint guns at different pressures. The optimal 
paint thickness was reported to be 1.0 mil for the inside of the grill and we consider it 
unacceptable if the average thickness exceeds that. We consider two factors, namely 
gun position and pressure, with two and three levels, respectively (for confidentiality, 
the actual levels were not reported). From Table 6.7 of Wu and Hamada (2009) the 
resulting experiment can be considered as a uniform design with six replicates and 
the data are given in Table 5.1. Consider a 2 x 3 main-effects model with logit link 
due to insignificant interactions. The design matrix is same as in expression ((2]). 
Fitting a logistic regression model, we get (3 = (0.234, -0.237, -0.341, -0.234)'. If 
one wants to conduct a follow-up experiment on half-fractions, then it is sensible to 
use the knowledge obtained by analyzing the data. Instead of taking (3 exactly if a 
point close to it is taken as the initial value, say, (3 = (0.25, —0.25, —0.30, —0.25)', the 
D-optimal design is p = (0.162, 0.178, 0.175, 0.175, 0.148, 0.162). 

Table 5.1: Paint Data (Hale-Bennett and Lin, 1997) 
Gun Pressure ^success out of 6 replicates 



1 1 2 

1 2 4 

1 3 3 

2 1 3 
2 2 4 
2 3 4 



For exact design, again we consider f3 = (0.25, —0.25, —0.30, —0.25)' as the true pa- 
rameter. Using the exchange algorithm with 36 experimental units, we obtained the 
D-optimal exact design is n = (6,7,6,6,5,6), which is close the uniform allocation 
(6, 6, 6, 6, 6, 6). 

Example 5.2: Printed Circuit Board: Jeng, Joseph and Wu (2008) reported an 
experiment on printed circuit boards (PCBs), which has been modified to suit our 
purpose. This experiment was about inner layer manufacturing of PCBs. Several 
faults can occur during the manufacturing of PCBs, of which shorts and opens in 
the circuit are the major ones. In this case, we consider the number of opens as 
response and two factors at two and three levels respectively. The factors are whether 
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preheating was done or not, and lamination temperature (95°C, 105°C and 115°C). 
Table 5.2 is obtained from Table 2 of Jeng, Joseph and Wu (2008). 



Table 5.2: PCB Data (Jeng, Joseph and Wu, 2008) 



Preheat 


Temperature 


^success out of 480 replicates 


1 


1 


120 


1 


2 


16 


1 


3 


25 


2 


1 


50 


2 


2 


51 


2 


3 


22 



Again, as before, the design matrix is given by expression fl2]). Fitting a logistic 
regression model, we get (3 = (—2.37,0.154,0.717,0.113)'. Let us assume that the 
true (3 as (-2.5,0.15,0.70,0.10)', the D-optimal design is p D = (0.216,0.186, 0.198, 
0.206, 0.115, 0.080). Applying the exchange algorithm for integer- valued allocations, 
we get the optimal allocations as (621, 535, 569, 593, 331, 231)' which is far away from 
being uniform, each with 480 replicates. 

6 Discussion 

In this paper, we extended the results of Yang, Mandal and Majumdar (2012b) to 
more general cases under generalized linear model setup. Our results allow the re- 
sponse could belong to a single-parameter exponential family, which includes Bino- 
mial, Poisson, exponential distributions as special cases. Our results also allow a fairly 
arbitrary set of design points, which could come from combinations of multiple-level 
factors, or grid points of continuous covariates. 

Yang, Mandal and Majumdar (2012b) also proposed EW D-optimal designs which is 
much easier to compute than the usual Bayesian D-optimal designs and more robust 
than uniform designs. The same concept and techniques can also be extended to a 
single-parameter exponential family response with arbitrary but fixed and distinct 
design points. Their discussion on fractional factorial designs can also extended here. 
This could be a topic of future research. 
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